Mathematical Modeling and Computer Simulations

Laboratories

1 Analysis of the romance model of romance using Romeo and Juliet as an example.

Strogatz derived a simple model, based on a system of ordinary differential equations, to describe the dynamics of a romance between two people. We will consider a discrete version of the model using the Romeo and Juliet romance as an example. A deterministic, dynamic, linear model with discrete time looks as follows.

\[\begin{align*}R_{n+1} = a_{R}R_{n} + p_{R}L_{n}\\ L_{n+1} = a_{L}L_{n} + p_{L}R_{n} \end{align*}\]

Interpretation of coefficients:

We are simulating this model on a time interval of 100 days (\(T = [0,100], \Delta t = 1\)), so we will start by creating functions for the model and graphs.

First, the simulation will be carried out for the parameters: \[R_0 = L_0 = 1, a_R = 0.5, a_L = 0.7, p_R = 0.7, p_L = 0.9\]

We start with the assumption that both Romeo and Juliet are initially in a state of not very strong love for each other (\(R_0 = L_0 = 1\)). The parameters \(a_R\), \(a_L\) take values less than 1, therefore, if Romeo and Juliet were to respond only to the dynamics of their own feelings (model (1), case \(p_R = p_L = 0\)), then their feelings toward each other would become neutral. In our case, Romeo and Juliet respond to each other, but rather tentatively (\(p_R\) and \(p_L\) are positive, with the strength of Juliet’s reaction to Romeo’s feelings being significant over time). We see that with this combination, their affection increases exponentially and goes to infinity at \(n \rightarrow \infty\).

Simulation two is carried out for the parameters: \[R_0 = L_0 = 1, a_R = 1.0, a_L = 1.0, p_R = 0.2, p_L = -0.2\]

Here we also start with the assumption that both Romeo and Juliet are initially in a state of not very strong love for each other (\(R_0 = L_0 = 1\)). The parameters \(a_R\), \(a_L\) take values = 1 therefore, if Romeo and Juliet were to respond only to the dynamics of their own feelings (model (1), case \(p_R = p_L = 0\)), then their feelings toward each other would be constantly at = 1. In our case, Romeo and Juliet react to each other with the same strength in a slight way the result of which we see a mood swing, because when Romeo’s feelings for Juliet increase then hers weaken and vice versa. Here we have a divergence in time.

Second simulation for parameters: \[R_0 = L_0 = 1, a_R = 0.5, a_L = 0.8, p_R = 0.2, p_L = 0.5\]

Again, we start with the same assumption that both Romeo and Juliet are initially in a state of not very strong love for each other (\(R_0 = L_0 = 1\)). From the values of \(a_R, p_R\), we can read that Romeo’s love for Juliet will weaken under the assumption of positive initial love between Romeo and Juliet until at \(n \rightarrow \infty\) it approaches zero. In turn, from the parameters \(a_L, p_L\) we will get that at a positive value of Romeo and Juliet’s affection (initial) Juliet’s love will increase. In such a scenario, over time, Romeo’s affection will become neutral, and Juliet’s will not increase much.

2 Determination of equilibrium points and stability analysis using the Romeo and Juliet model as an example

## Ładowanie wymaganego pakietu: MASS

We will now try to understand the result of Romeo and Juliet romance observed in the four given cases from the point of view of stability of any fixed points. We consider again 3 cases of model parameters:

For simulation 1:

## [1] "Solution R*=0 and L*=0"
## [1] "Checking stability..."
## [1] "The determinant of a homogeneous system: -0.48"
## [1] "A homogeneous system has exactly one solution!"
## [1] "Trace value: 1.2 and the determinant value: -0.28"
## [1] "Jury conditions not met, point not stable!"

For simulation 2..:

## [1] "Solution R*=0 and L*=0"
## [1] "Checking stability..."
## [1] "The determinant of a homogeneous system: 0.04"
## [1] "A homogeneous system has exactly one solution!"
## [1] "Trace value: 2 and the determinant value: 1.04"
## [1] "Jury conditions not met, point not stable!"

For simulation 3..:

## [1] "Solution R*=0 and L*=0"
## [1] "Checking stability..."
## [1] "The determinant of a homogeneous system: 0"
## [1] "The system has an infinite number of fixed points!"

We see that for no point are Jury’s conditions satisfied. Indeed, if we look at the solutions they do not converge to a single fixed point (\(R^{*}=0,L^{*}=0\)).

Simulation 3, on the other hand, contains infinitely many fixed points. We will determine, therefore, their family using the equation \[\begin{align*} (a_r - 1)\cdot R^{*} + p_R\cdot L^{*} = 0 \\ (a_L - 1)\cdot L^{*} + p_L\cdot R^{*} = 0 \end{align*}\] Because we can say that the total feeling of love/hate between Romeo and Juliet is initially equal to the \(R_0 + L_0\) and preserved throughout the observation period, then we can assume that the \(R^{*} + L^{*} = R_0 + L_0\). Because of the initial conditions \(R_0 = 1\) and \(L_0=1\) we have \(L^{*} = 2 - R^{*}\). Thus, taking the initial equation, we have a family of fixed points \[\{(R^{*}, L^{*}) : R^{*} = \frac{2\cdot p_R}{(1- a_R)+p_R}, L^{*} = 2-R^{*}\}\].

Thus, for this case, the exact value of the fixed point will be

## [1] "Wartość R*= 0.571428571428572"

## [1] "Wartość L*= 1.42857142857143"

We now turn to determining the eigenvalues of the Jacobian. We have a Jacobian equal to \[ J = \begin{bmatrix} 0.5 & 0.2 \\ 0.5 & 0.8 \end{bmatrix} \] The characteristic equation is thus of the form: \[ \lambda^2 - trJ \cdot \lambda + detJ = 0 \] Inserting the corresponding values \(trJ = 1.3\) and \(detJ = 0.3\) we have: \[\lambda^2 - 1.3\cdot\lambda + 0.3 = 0\] skąd \[(\lambda-1)(\lambda-0.3) = 0\] so \(\lambda_1 = 1\) and \(\lambda_2 = 0.3 = 0.5 + 0.8 - 1 = a_r + a_L - 1\) for this specific case. From a theorem, we obtain that the point we have determined is stable by virtue of \(|\lambda_2| < 1\). It is important to note that the possibility of such an inference with one eigenvalue satisfying the conditions is possible only for linear systems.

3 Prey-predator type models.

We will consider a discrete model of the Lotki-Volterra type assuming unlimited capacity of the environment for the prey population. Suppose we have populations of two species, e.g. foxes (predator - predator ) and rabbits (prey - prey). Let \(P_n\) denote the density of the rabbit population, and \(Q_n\) denote the density of the fox population at time \(n\). The predator-prey model has the following form: \[\begin{align*} P_{n+1} = P_n(1 + r(1 - \frac{P_n}{K})) - sP_nQ_n = f_1(P_n, Q_n) \\ Q_{n+1} = (1-u)Q_n + asP_nQ_n = f_2(P_n, Q_n) \end{align*}\] where

Assume that \(P_n = P_{n+1} = P^{*}\) and \(Q_n = Q_{n+1} = Q^{*}\) then the system \[\begin{align*} P^{*} = f_1(P^{*}, Q^{*})\\ Q^{*} = f_2(P^{*}, Q^{*}) \end{align*}\] has three solutions (equilibrium points) of the form \[\begin{align*} (P^{*}_1, Q^{*}_1) = (0,0) \\ (P^{*}_2, Q^{*}_2) = (K, 0) \\ (P^{*}_3, Q^{*}_3) = (\frac{u}{sa}, \frac{r}{s}(1-\frac{u}{saK})) \end{align*}\] Indeed, from the second equation of the system we get the following \[\begin{align*} Q^{*} = Q^{*}[(1-u) + asP^{*}] \\ Q^{*}[u - asP^{*}] = 0 \end{align*}\] from where \(Q^{*}=0\) or \(P^{*} = \frac{u}{as}\). Let us first consider the first situation \(Q^{*} = 0\), then by inserting into the first equation of the system we have: \[\begin{align*} P^{*} = P^{*}(1 + r(1 - \frac{P^{*}}{K})) \\ P^{*} = P^{*} + P^{*}r - \frac{(P^{*})^2r}{K} \\ P^{*}r - \frac{(P^{*})^2r}{K} = 0 \\ P^{*}(r - \frac{P^{*}r}{K}) = 0 \end{align*}\] so \(P^{*} = 0\) or \(P^{*} = K\). Thus, we have already demonstrated the equality of existence of two equilibrium points \[\begin{align*} (P^{*}_1, Q^{*}_1) = (0,0) \\ (P^{*}_2, Q^{*}_2) = (K, 0) \end{align*}\]. Now we will consider the situation where \(P^{*} = \frac{u}{as}\). Inserting into the first equation of the system we have: \[\begin{align*} \frac{u}{as} = \frac{u}{as} + \frac{ur}{as} - \frac{u^2r}{(as)^2K} - \frac{u}{a}Q^{*} \\ Q^{*} = \frac{r}{s} - \frac{ur}{as^2K} \\ Q^{*} = \frac{r}{s}(1 - \frac{u}{saK}) \end{align*}\] From there we get the third equilibrium point \[ (P^{*}_3, Q^{*}_3) = (\frac{u}{sa}, \frac{r}{s}(1-\frac{u}{saK})). \]

We create a function for the model and graphs.

We will now simulate the model for the following initial parameters: \[ N = 100, P_0 = 0.1, Q_0 = 0.2, r=0.5, u = 0.2, s = 2.2, a = 0.9, K = 0.7 \]

Let’s check the stability of the three equilibrium points for the mapping \(f=(f_1, f_2)\). The Jacobi matrix then has the form: \[ J = \begin{bmatrix} 1 + r - \frac{2r}{K}P^{*} - sQ^{*} & -sP^{*} \\ asQ^{*} & 1 - u + asP^{*} \end{bmatrix} \]

For the parameters and individual equilibrium points given above, we get:

From the simulation for \(N=100\) is the stability of this point is seen in a rather weak way. Moreover, it looks like cyclic stability. So let’s check the simulation for \(N=1000\) generation.

Here we can already see the stability of the indicated point with its cyclicality. Moreover, from the phase portrait we can clearly see the pursuit of the indicated stable point.

In the given simulation, the initial density of rabbits is twice as low as the initial density of foxes. What’s more, foxes have a high attack rate, so in the first generations we can see a large decrease in rabbit population (higher fox population and high effectiveness). However, after the rabbit population declines, the predators run out of food and their population numbers also decline. Meanwhile, the population of prey (rabbits) begins to increase again, and the foxes have food - another attack occurs. During the attack, the population of foxes increases and the population of rabbits decreases. This turn of events explains the cyclical nature of the indicated equilibrium point.

4 Continuous models and the problem of choosing an appropriate numerical method

In what follows, we consider the predator-prey model in the context of two differential equations (continuous type model), i.e. \[\begin{align*} \frac{dP}{dt} = rP(1-\frac{P}{K}) - sPQ = f_1(P, Q) \\ \frac{dQ}{dt} = -uQ + asPQ = f_2(P,Q) \end{align*}\] The discrete model considered in the previous lab is a unit time map for the above system of differential equations. The equilibrium points of this system do not change therefore \[ (P^{*}, Q^{*}) = (0,0),\quad (P^{*}, Q^{*}) = (K,0),\quad (P^{*}, Q^{*}) = (\frac{u}{sa}, \frac{r}{s}(1-\frac{u}{saK})) \] The Euler method for this system has the form \[\begin{align*} P_{n+1} = P_n + \Delta t(rP_n(1-\frac{P_n}{K}) - sP_nQ_n), \\ Q_{n+1} = Q_n + \delta t(-uQ_n + asP_nQ_n), \end{align*}\] Where at \(\Delta t = 1\) we get the discrete system from the previous sections.

In turn, the RK2 method for this system takes the form:

\[\begin{align*} k_{1P}&= f_1(P_n, Q_n)\\ k_{1Q}&= f_2(P_n, Q_n)\\ k_{2P}&= f_1(P_n + \Delta t k_{1P}, Q_n + \Delta t k_{1Q}) \\ k_{2Q}&= f_2(P_n + \Delta t k_{1P}, Q_n + \Delta t k_{1Q}) \end{align*}\]

\[ \begin{cases} P_{n+1} = P_n + \frac{\Delta t}{2}(k_{1P} + k_{2P}),\\ Q_{n+1} = Q_n + \frac{\delta t}{2}(k_{1Q} + k_{2Q}) \end{cases} \]

The RK4 method for this system is of the form:

\[\begin{align*} k_{1P}&= f_1(P_n, Q_n) \\ k_{1Q}&= f_2(P_n, Q_n) \\ k_{2P}&= f_1(P_n + \frac{\Delta t}{2}k_{1P}, Q_n + \frac{\Delta t}{2}k_{1Q})\\ k_{2Q}&= f_2(P_n + \frac{\Delta t}{2}k_{1P}, Q_n + \frac{\Delta t}{2}k_{1Q})\\ k_{3P}&= f_1(P_n + \frac{\Delta t}{2}k_{2P}, Q_n + \frac{\Delta t}{2}k_{2Q})\\ k_{3Q}&= f_2(P_n + \frac{\Delta t}{2}k_{2P}, Q_n + \frac{\Delta t}{2}k_{2Q})\\ k_{4P}&= f_1(P_n + \frac{\Delta t}{2}k_{3P}, Q_n + \frac{\Delta t}{2}k_{3Q})\\ k_{4Q}&= f_2(P_n + \frac{\Delta t}{2}k_{3P}, Q_n + \frac{\Delta t}{2}k_{3Q}) \end{align*}\]

\[ \begin{cases} P_{n+1} = P_n + \frac{\Delta t}{6}(k_{1P} + 2k_{2P} + 3k_{3P} + 4k_{4P}), \\ Q_{n+1} = Q_n + \frac{\delta t}{6}(k_{1Q} + 2k_{2Q} + 3k_{3Q} + 4k_{4Q}). \end{cases} \]

We assume the following initial parameters: \[P_0 = 0.3, Q_0 = 0.9, r=1.2, u = 0.2, s = 2, a = 0.5, K = 1. \]

We will examine the stability of the equilibrium points for these points. For the parameters and individual equilibrium points given above, we will obtain:

For the indicated initial parameters, we will perform the following simulations (\(N=\frac{T}{\Delta t}\))and present phase portraits of the solutions. First, for the following values:

1. \(T = 100, \Delta t=0.1\).

According to the above simulations, for a time step of \(\Delta t = 0.1\), we can observe the convergence of the solutions of all methods to the equilibrium point of the system \((0.2,0.48)\). The methods give very close solutions, and RK2 and RK4 even coincide. Of course, with the numerical parameters chosen in this way, the assumptions about their convergence are met for each method, which are strictly dependent on the value of the time step. It can be seen that we obtain approximate solutions already on the interval \(t\in [0,30]\). We will now extend the number of days and increase the time step.

2. \(T = 1000, \Delta t=1\).

The above simulations for a time step of \(\Delta t = 1\) show the convergence of the solutions of the RK2 and RK4 methods to the equilibrium point of the system \((0.2,0.48)\), while although the Euler method also converges to this point, it falls into oscillations by which it pursues the solution in a less stable manner and requires more iterations. Still, the RK2 and RK4 methods even overlap after a small amount of time. Since the Euler method requires a small time step \(\Delta t\) let’s increase it slightly and check the results.

3. \(T = 1000, \Delta t=1.1\).

We can see that for a time step of \(\Delta t = 1.1\) therefore an increase from the previous simulation of just \(0.1\) the Euler method is already diverging, while the RK2 and RK4 methods continue to converge to the equilibrium point of the system \((0.2,0.48)\). This allows us to see the aforementioned instability of Euler’s method and a sizable dependence on the choice of time step. Let’s increase the time step again, this time by three times.

4. \(T = 1000, \Delta t=3\).

For a time step of \(\Delta t = 3\), only the RK4 method already converges to the equilibrium point of the \((0.2, 0.48)\) system, where both the Euler and RK2 methods already diverge. Moreover, both methods violate the basic assumption of non-negativity of population density. In this case, the RK4 method is the most robust to changing the time step size. Let’s test its further capabilities and again increase the step by a small amount.

5. \(T = 1000, \Delta t=3.2\).

In such a variant, no method converges to the equilibrium point of the solution of a given system anymore. Since we do not always have an influence on the choice of the simulation step because, for example, it can be imposed by the phenomenon itself, this shows the difficulty of modeling and choosing the right method to determine the solutions of certain equations. Each of them has assumptions to which special attention should be paid. As can be seen, even the assumption of a small time step size can significantly disturb the results and not lead to convergence of the solution.

5 Custom Mickens-type numerical model.

In what follows, we consider the predator-prey type model in the context of two differential equations (continuous type model), i.e.

\[\begin{align*} \frac{dP}{dt}&= rP(1-\frac{P}{K}) - sPQ = f_1(P, Q) \\ \frac{dQ}{dt}&= -uQ + asPQ = f_2(P,Q) \end{align*}\]

Such a mathematical model described by a differential problem, we are able to write in the form: \[\frac{dx}{dt} = f_+(x) - xf_-(x)\] where \(f_+, f_-: \mathbb{R}^m_{+} \bigcup \{0\} \rightarrow \mathbb{R}^m_{+} \bigcup \{0\}\) are \(C^1\) class. We can approximate such a model with a non-standard differential scheme of the form: \[ x_{n+1} = \frac{x_n + \phi(\Delta t; p)f_+(x_n)}{1 + \phi(\Delta t; p)f_-(x_n)} \] Such a numerical model is unconditionally convergent regardless of the time step \(\Delta t\) unlike previous numerical methods as shown in labs 4 and 5. When obtaining functions \(f_+\) and \(f_-\) with positive values and starting from a non-negative value \(x(t_0) = x_0 \geq 0\) guarantee non-negative values. We can write our dynamic system as

\[\begin{align*} \frac{dP}{dt}&= rP - P(r\frac{P}{K} + sQ)\\ \frac{dQ}{dt}&= asPQ - Qu \end{align*}\]

Note that we can write the first equation in the desired form: \[\frac{P_{n+1} - P_n}{\Delta t} = rP_n - P_{n+1}(r\frac{P_n}{K} + sQ_n).\] We can transform such an equation to the following discrete system assuming a trivial step function \(\phi(\Delta t;p) = \Delta t\): \[P_{n+1} = \frac{P_n + \Delta t\cdot rP_n}{1 + \Delta t\cdot (r\frac{P_n}{K} + sQ_n)}\] On the other hand, we can write the second equation in the form \[\frac{Q_{n+1} - Q_n}{\Delta t} = asP_nQ_n - Q_{n+1}u\] and transforming as: \[Q_{n+1} = \frac{Q_n + \Delta t\cdot (asP_nQ_n)}{1 + \Delta t\cdot (u)}\] Let’s implement these numerical methods and check the solutions.

In the first place for the parameters:

1. \(T = 1000, \Delta t=0.1\).

We see that for the given time step the discretization of the initial dynamical system is correct. Simulation-wise, we get the desired solutions, i.e. P and Q are moving towards a stable equilibrium point \((0.2,0.48)\). Let’s check the larger time step.

2. \(T = 1000, \Delta t=1.1\).

We can see that although the solutions strive for a stable equilibrium point they fall into oscillatory motion and are not stable. It takes a large number of iterations for them to reach convergence. Let’s see what happens for an even larger time step, assume:

3. \(T = 1000, \Delta t=3.2\).

We see that now the solutions have not reached convergence at all despite the discretization. Looking for the reason for this behavior, we will come to the way the discretization of the given dynamical system was implemented. It is incorrect because it is two processes that have been discretized separately and do not interact with each other. So we will make a minor change in the program so that the process \(Q\) is dependent at time \(n\) on the process \(P\) also at time \(n\), rather than at \(n-1\). So we assume the following dynamic system after discretization: \[P_{n+1} = \frac{P_n + \Delta t\cdot rP_n}{1 + \Delta t\cdot (r\frac{P_n}{K} + sQ_n)}\] \[Q_{n+1} = \frac{Q_n + \Delta t\cdot (asP_{n+1}Q_n)}{1 + \Delta t\cdot (u)}\]

Let’s check again for the time step

2. \(T = 1000, \Delta t=1.1\).

Here we see that convergence is achieved much faster than before. Let’s see if it will be achieved at all for the next example.

3. \(T = 1000, \Delta t=3.2\).

As you can see, the solutions reach convergence to the desired equilibrium point. The second implementation includes the compatibility of the dynamics of the two processes according to continuous time therefore reaches convergence to the solutions. Let’s see for an even larger time step.

4. \(T = 1000, \Delta t=5\).

Convergence has been achieved and in this case, but it can be seen that some small oscillations of the solutions are introduced so more iterations are probably already needed. Let’s assume, then, some arbitrarily large step and test the method.

5. \(T = 1000, \Delta t=10\).

For such a large time step, we can see the slowness of reaching the stability point. Convergence is reached after more than 1000 iterations. Nevertheless, the discretization still maintains stability and the system converges to a solution.

Let’s also try to find a non-trivial time step function \(\phi(\Delta t;p)\) that will help describe the dynamics in our equation. To do this, let’s return to our initial equation \[ \frac{dP}{dt} = rP - P(r\frac{P}{K} + sQ) = rP(1 - \frac{P}{K}) + sPQ \] Note that it is similar to a logistic equation of the form \[\frac{dy}{dt} = \lambda y(1-y)\] where we would only have a simple scaling relative to our equation. Thus, we can expect similar dynamics of the two equations and use the construction of a non-trivial step function for the logistic equation. In this case, we consider (using the implicit Euler method) a numerical solution of the form \[y_{k+1} = \frac{1}{1 + \lambda h}y_k,\] where the significant factor depending on the parameter \(1+ \lambda h\) can be written in the form \[1 + \lambda h = e^{\lambda h} + O(\lambda^2 h^2)\] and with a sufficiently small \(0 < \lambda h \ll 1\) is assumed. \[ h = \frac{e^{\lambda h} - 1}{\lambda} \] so the non-trivial step function is of the form \(\phi(h) = \frac{e^{\lambda h} - 1}{\lambda}\). In our situation, \(h = \Delta t\) and \(\lambda = r\), so the non-trivial step function for the first equation from the dynamical system is the function \[ \phi(\Delta t, r) = \frac{e^{r\Delta t} - 1}{r} \] and inserting into the exact differential diagram we get \[P_{n+1} = \frac{P_n + (\frac{e^{r\Delta t} - 1}{r})rP_n}{1 + (\frac{e^{r\Delta t} - 1}{r})(r\frac{P_n}{K} + sQ)}.\] Let’s now perform analogous reasoning for the second equation of our dynamic system. We have \[ \frac{dQ}{dt} = -uQ + asPQ = Q(asP-u), \] and so again a scaled version of the logistic equation. Thus, following the same steps, we obtain a non-trivial step function where \(\lambda = u\) of the form: \[\phi(\Delta t, u) = \frac{e^{u\Delta t} - 1}{u}\] and inserting into the exact differential diagram we get

\[Q_{n+1} = \frac{Q_n + (\frac{e^{u\Delta t} - 1}{u})(asP_nQ_n)}{1 + (\frac{e^{u\Delta t} - 1}{u})u}\] Let’s see the results using the same examples.

We can see that using a non-trivial step function convergence is also achieved, but not at all in a faster way. For a large time step, a large number of iterations is further needed. This shows the importance of the correctness of the implementation of the dynamic model to obtain convergence of solutions for different time steps.

Comparing to Euler’s method, we can observe, first of all, the stability of the solutions, since Euler’s method ceased to be stable already for a time step of \(\Delta = 1\). On the other hand, the alternative methods RK2 and RK4, although they were more stable (however, for a larger time step they also stopped converging), in addition to this, in the case of a multi-equation model (we only have 2) the implementation would be extremely difficult due to the need to write many equations. This argues for a discretization of the dynamical system.

The Jacobian matrix for our dynamical system is of the form: \[ J = \begin{bmatrix} 1 + r - \frac{2r}{K}P^{*} - sQ^{*} & -sP^{*} \\ asQ^{*} & 1 - u + asP^{*} \end{bmatrix} \] and for the equilibrium point \((P^{*}, Q^{*}) = (0,0)\) it is as follows: \[ J = \begin{bmatrix} 1 + r & 0 \\ 0 & 1 - u \end{bmatrix} \]

Assuming the parameters from the task, we get \[ J= \begin{bmatrix} 2.2 & 0 \\ 0 & 0.8 \end{bmatrix} \] Then we get \(trJ = 3\) and \(detJ = 1.76\) so Jury’s conditions are not satisfied because \(trJ = 3 > 2.76 = 1 + detJ\) and \(1+detJ = 2.76 > 2\). Thus, the point \((0,0)\) is not stable for the parameters we have assumed. This is logical, because with an initial population of rabbits and foxes (prey and predators) with some density and positive growth rates for the prey population and fox birth rate, we do not expect both populations to die out (with these parameters and this model design). Of course, these populations in the real world may go extinct with some random factor, but in our situation the model is not able to predict such a situation.

We will now check for the given dynamical system the built-in deSolve package that allows solving initial problems for differential equations/systems.

Let’s first check the solutions with our parameters for the lsoda function.

## 
## Dołączanie pakietu: 'deSolve'

## Następujące obiekty zostały zakryte _przez_ '.GlobalEnv':
## 
##     euler, rk4

We see indeed a fast convergence to the equilibrium point and stability of solutions. Let’s check for the ode function.

And here, too, we see seamless convergence to the equilibrium point of the embedded function. Let us now juxtapose the solutions of these methods with the methods tested by us.

6 Parameter estimation and population growth model

6.1 Parameter estimation

6.1.1 Theoretical part

We assume that \(X_1, ..., X_n\) are independent random variables and the observed data \(x_i\), \(i=1,...,n\) are realizations of random variables \(X_i\) with Pareto distribution with density function \[ f(x_i | x_0, \lambda) = \lambda x_0^\lambda x_i^{-\lambda-1}, \] Where \(x_i \geq x_0\), > 1$, \(x_0\) fixed. Then the credence function for the parameter \(\lambda\) has the following form (since these variables are independent then the product of the density function is simply their product): \[ \mathcal{L}(\lambda) = \prod_{i=1}^{n}{f(x_i|x_0, \lambda)} = \prod_{i=1}^{n}{\lambda x_0^\lambda x_i^{-\lambda-1}} = \lambda^n x_0^{n\lambda}\prod_{i=1}^{n}{x_i^{-\lambda-1}} = \lambda^n x_0^{n\lambda}(x_1^{-\lambda - 1}\cdot ... \cdot x_n^{-\lambda - 1}). \]

Taking the logarithm of the credibility function, we get: \[ \mathcal{L}\mathcal{L}(\lambda) = \ln \mathcal{L}(\lambda) = \ln (\lambda^n x_0^{n\lambda} \prod_{i=1}^{n}{x_i^{-\lambda -1}}) = \ln (\lambda^n x_0^{n\lambda}(x_1^{-\lambda-1}\cdot ... \cdot x_n^{-\lambda-1})) \] We are looking for the estimate of the highest reliability of the parameter \(\lambda\) so we need to find the point \(\overline{\lambda}\) for which the function \(\mathcal{L}\mathcal{L}(\lambda)\) reaches its extremum. Thus, we are looking for the solution of the equation \(\frac{d\mathcal{L}\mathcal{L}(\lambda)}{d\lambda} = 0\). Let’s start by transforming the function \(\mathcal{L}\mathcal{L}(\lambda)\) i.e.: \[ \mathcal{L}\mathcal{L}(\lambda) = \ln (\lambda^n x_0^{n\lambda}(x_1^{-\lambda-1}\cdot ... \cdot x_n^{-\lambda-1})) = \ln{\lambda^n} + \ln{x_0^{n\lambda}} + \sum_{i=1}^{n}{\ln{x_i^{-\lambda-1}}} = n\ln{\lambda} + n\lambda\ln{x_0} + (-\lambda - 1)\sum_{i=1}^{n}{\ln{x_i}} \] Now taking the derivative after \(\lambda\) we get: \[ \frac{d\mathcal{L}\mathcal{L}(\lambda)}{d\lambda} = \frac{d}{d\lambda}(n\ln{\lambda} + n\lambda\ln{x_0} + (-\lambda - 1)\sum_{i=1}^{n}{\ln{x_i}}) = \frac{d}{d\lambda}n\ln{\lambda} + \frac{d}{d\lambda}n\lambda\ln{x_0} + \frac{d}{d\lambda}(-\lambda -1)\sum_{i=1}^{n}{\ln{x_i}} = \frac{n}{\lambda} + n\ln{x_0} - \sum_{i=1}^{n}{\ln{x_i}} \] and comparing \(\frac{d\mathcal{L}\mathcal{L}(\lambda)}{d\lambda} = 0\) we get: \[ \frac{n}{\lambda} + n\ln{x_0} - \sum_{i=1}^{n}{\ln{x_i}} = 0 \] \[ \frac{n}{\lambda} = \sum_{i=1}^{n}{\ln{x_i}} - n\ln{x_0} \] \[ \lambda = \frac{1}{\frac{1}{n}\sum_{i=1}^{n}{\ln{x_i}} - \ln{x_0}} \] Thus, we obtain that the estimate of the highest reliability of the parameter is \[ \overline{\lambda} = \frac{1}{\frac{1}{n}\sum_{i=1}^{n}{\ln{x_i}} - \ln{x_0}}. \]

As initially thought.

6.1.2 Simulation part

We will proceed to the simulation part. To do this, we will load the file , which contains a 20-element sample with a Pareto distribution generated with a density function with parameter \(\lambda = 1.5\). We will also define a function \(-\mathcal{L}\mathcal{L}(\lambda)\), which we will minimize (equal to maximizing the function \(\mathcal{L}\mathcal{L}\)).

Let’s now plot the defined function for our sample and parameters \(x_0 = 0.9\) and \(\lambda \in [1,3]\).

Using the function we will minimize the function \(-\mathcal{L}\mathcal{L}\) in order to obtain an estimate of the optimal parameter \(\overline{\lambda}\).

l_0 <- 1.1
min_l <- nlm(f=ll, p= l_0, x_0=x_0, x =x)

Thus, we obtained an estimate of the minimum of the function equal to \(21.44467\) for the parameter \(\overline{\lambda} = 1.577962\). It is close to the known parameter \({\lambda} = 1.5\) of the Pareto distribution of the sample. We will now compare it with the theoretical estimate by us derived for the same data sample.

n <- length(x)

min_l_t <- 1 / (1/n * (sum(log(x))) - log(x_0))
min_ll <- ll(min_l_t, x_0, x)

In this case, we have a minimum value equal to \(21.44467\) for the parameter \(\overline{\lambda} = 1.577962\), so we see a correspondence between the simulation solution and the theoretical value.

Now let’s check the results for the initial values \((x_0, \lambda_0) = (0.85, 1)\).

l_0 <- 1
x_0 <- 0.85

min_l <- nlm(f=ll, p= l_0, x_0=x_0, x =x)

We get the minimum value of the function equal to \(23.17178\) and an estimate of the parameter \(\overline{\lambda}=1.447414\), which is still close to the actual value of the parameter. Let’s check the theoretical estimation for the new starting point.

min_l_t <- 1 / (1/n * (sum(log(x))) - log(x_0))
min_ll <- ll(min_l_t, x_0, x)

We get a minimum value of \(23.17178\) and an estimate of the parameter \(\overline{\lambda}=1.447415\), so again the results are the same. Consider the other end of the \(\lambda\) interval along with a new starting point \((x_0,\lambda_0) = (0.95, 3)\).

l_0 <- 3
x_0 <- 0.95

min_l <- nlm(f=ll, p= l_0, x_0=x_0, x =x)

Now the function based on the Newton-type optimization algorithm has not reached convergence. The minimum value of the function determined is equal to \(-709.7826\), and the estimate is \(-55.21543\). Let’s see the results for the theoretical solution.

min_l_t <- 1 / (1/n * (sum(log(x))) - log(x_0))
min_ll <- ll(min_l_t, x_0, x)

We get the value of the function equal to \(19.66114\) with the estimation of the parameter \(\overline{\lambda} = 1.725145\). While the estimation deviates more from the actual value of the parameter, the error is still negligible. The Newton-type optimization algorithm is based on the expansion of the function into a second-order Taylor series hence the possibility of divergence when the initial values are incorrectly chosen. We will simulate the potential parameter estimates depending on the choice of the initial \(\lambda_0\) value to see what effect this has on the optimization result.

As speculated, the initial value of \(\lambda_0\) affects the obtained optimization results. As the value increases, the algorithm becomes divergent. Constraints should be imposed on the optimization of the parameter, which are available in other packages. This simulation shows the importance of the choice of initial parameters.

6.2 Population growth model

Let us denote by \(N_D\) the vector of observed population abundances for the specified years. By \(N(t)\) we will denote the population size at time \(t \geq 0\). The dynamics of population change can be described by Verhulst’s logistic growth equation, we assume the model: \[ \frac{dN(t)}{dt} = rN(t)(1-\frac{N(t)}{K}), \quad N(0)=N_0 \] Where \(N_0 > 0\) is the initial value of the population size, \(r>0\) the population growth parameter and \(K>0\) is the capacity of the environment. Such a vector of parameters takes the form \(p=\{N_0, r, K\}\). We determine the initial condition \(N_0\) from the data. Solving the above problem, we get a solution of the form \[ N(t, p) = \frac{N_0 K}{N_0 + (K-N_0)e^{-rt}}. \] We want to minimize the error of least squares estimation of population size, which is of the following form \[ BNK(p) = \sum_{i=1}^{n}{(\overline{N_D} - N(t_i, p))^2}, \] where \(\overline{N_D} = \frac{N_D}{1000}\) is the scaled value (you can check the behavior without scaling). In , we can use the functions to scale and the function to minimize the BNK from where we get the parameter estimate.

We load a file containing the population size of Gdansk from 1950 to 2015 with a 5-year interval. We estimate by BNK minimization the \(p\) parameters of the population model based on this data assuming \[ \overline{N_D} = \frac{N_D}{10000}, \] \[ p = (\overline{N_0}, r, K) = (\overline{N_D}(1), 0.01, 2\cdot \max{(\overline{N_D})}). \]

From the minimization we obtained the optimal parameters equal to \[ \overline{N_0} = 19.463209, \] \[ \overline{r} = -2.757316, \] \[ \overline{K} = -1424.408053. \]

The solution fit on these estimated parameters obtained using least squares minimization is as follows:

However, let’s note that despite the fact that the estimated parameters do not meet the assumptions, the model fit is quite good. Now let’s check the behavior of the algorithm for a smaller scaling with the same initial parameters ie. \[ \overline{N_D} = \frac{N_D}{100}, \] \[ p = (\overline{N_0}, r, K) = (\overline{N_D}(1), 0.01, 2\cdot \max{(\overline{N_D})}). \]

From the minimization we obtained the optimal parameters equal to \[ \overline{N_0} = 1946.33, \] \[ \overline{r} = 0.1236176, \] \[ \overline{K} = 9372.32 \]

We can see that while the estimated value of the initial population value has been the same only scaled, the other parameters have changed. The solution fit on these estimated parameters obtained using least squares minimization looks as follows:

As can be seen, the model fit with the new estimated parameters with a smaller scale is much worse. Evidently, this is affected by the overestimated values of the estimated parameters \(\overline{r}\) and \(\overline{K}\), which are much larger and correspond to the scale of the population size. However, the estimated parameters meet the model’s assumptions of positivity. So let’s check again for a larger scale with the same parameters ie. \[ \overline{N_D} = \frac{N_D}{1000000}, \] \[ p = (\overline{N_0}, r, K) = (\overline{N_D}(1), 0.01, 2\cdot \max{(\overline{N_D})}). \]

From the minimization we obtained the optimal parameters equal to \[ \overline{N_0} = 0.1338027, \] \[ \overline{r} = 0.4714063, \] \[ \overline{K} = 0.4734332 \]

Now, in turn, the estimation of the initial population value parameter has decreased, and the other parameters have increased. The value of the \(\overline{K}\) parameter is close to the value of the \(\overline{r}\) parameter. The solution fit on these estimated parameters obtained using least squares minimization is as follows: The model fit for such an estimation is similar to that for the initial scale. Moreover, the estimated parameters are positive so they are consistent with the model assumptions. This indicates that the Newton-type optimization algorithm is more sensitive for functions taking high values and it is more difficult to obtain optimal parameter estimates in such situations. This algorithm is based on iterative counting of argument values using the Taylor series of a function of order two, so naturally, at high values of this function, the estimates may contain a large error threshold.

So let’s still try for the same scaling to take other starting parameters, i.e.  \[ \overline{N_D} = \frac{N_D}{10000}, \] \[ p = (\overline{N_0}, r, K) = (\overline{N_D}(1), 1, 10\cdot \max{(\overline{N_D})}). \]

From the minimization we obtained the optimal parameters equal to \[ \overline{N_0} = -102.9745 \] \[ \overline{r} = -80.09881 \] \[ \overline{K} = -78.44528 \]

Our parameters are inconsistent with the assumptions. Even the initial value of the population size is negative. Let’s check the situation if we reduce the initial parameters accordingly ie: \[ \overline{N_D} = \frac{N_D}{10000}, \] \[ p = (\overline{N_0}, r, K) = (\overline{N_D}(1), 0.0001, 0.1\cdot \max{(\overline{N_D})}). \]

From the minimization we obtained the optimal parameters equal to \[ \overline{N_0} = 0.2972155 \] \[ \overline{r} = -0.001612179 \] \[ \overline{K} = 0.1506714 \]

So this time further the parameter \(\overline{r}\) is inconsistent with the assumption. Nevertheless, let’s plot the graph of the function with these parameters.

It can be seen that it is indeed not fitted to the data and the results are not accurate. It turns out, then, that the best combination is to adopt high scaling and properly adjusted initial parameters. Scaling too low can lead to incorrect parameter estimates that are inconsistent with the assumptions, as can adopting inappropriate initial conditions.

7 Simulation of Brownian motion

7.1 Generating the trajectory of the process

A time-dependent function whose values at each time step are random variables is called a stochastic process.

A stochastic process of the form

\[ \{B(t), \quad t\in [0,T]\} \]

On the other hand, we call a Brownian motion process (Wiener process) if it satisfies the following conditions - the increments of the process are independent random variables, i.e. for any moments \(t_0 < t_1 < ... < t_n \leq T\) increments.

\[ B(t_1) - B(t_0), B(t_2) - B(t_1), ..., B(t_n) - B(t_{n-1}) \]

are independent; - the process increments are homogeneous (stationary), i.e:

\[ \forall_{t>s} \quad B(t) - B(s) \sim \mathcal{N}(0, t-s); \]

• The process has continuous trajectories with probability 1 ie:

\[ P(\{\omega \in \Omega : t \mapsto B(t,\omega) \textrm{ jest funkcją ciągłą}\}) = 1; \]

If in addition

\[ P(B(0) = 0) = 1 \]

then the stochastic process \(\{B(t), t \in [0, T]\}\) is called a standard Brownian motion process. We consider a discretized Brownian process i.e. one defined on a finite, predetermined set of division points of the segment \([0,T]\). Let us note that based on the properties of Brownian motion for a fixed increment \(\Delta t >0\), it occurs

\[ B(\Delta t) = B(\Delta t) - B(0) \sim \mathcal{N}(0, \Delta t) \sim \sqrt{\Delta t}\mathcal{N}(0, 1) \]

and

\[ B(t + \Delta t) - B(t) \sim \mathcal{N}(0, \Delta t) \sim \sqrt{\Delta t}\mathcal{N}(0, 1). \]

7.2 Algorithm for generating trajectories.

For the interval \([0, T]\) we set \(N>0\) and take the time step \(\Delta t = \frac{T}{N}\) obtaining the following segmentation of the segment.

\[ 0 = t_0 < t_1 < ... < t_{N-1} < t_N = T \quad t_i = i\Delta t. \]

Let us denote by \(B_i = B(t_i)\) the value of the trajectory of the Brownian motion process at the point \(t_i\). Then we generate the trajectory of the process according to the following algorithm: 1. \(B_0 = 0\); 2. \(B_i = B_{i-1} + \sqrt{\Delta t}\epsilon_i\), where \({\epsilon_i, \quad i=1, ..., N}\) is a sequence of independent random variables with distribution \(\mathcal{N}(0,1)\). Such a simulation generates trajectory points only at the dividing points of the segment \([0, T]\) between the dividing points the trajectory is approximated by linear interpolation (in the program we connect the determined discrete points with lines).

An example function that generates a trajectory of Brownian motion depending on the parameters \(N\) and \(T\) looks as follows:

brown <- function(N, T){
  return(cumsum(c(0, sqrt(T/N) * rnorm(N, 0, 1))))
}

Now let’s generate process trajectories for \(T=1\) with three different segment divisions, i.e., we will assume \(N \in \{2^5, 2^7, 2^9\}\). The results are as follows

We can see how much influence a change in the value of \(N\) has on the number of “bounces” of the Brownian motion trajectory. Moreover, a small number of bounces that is, at the same time, a high number of \(\Delta t\) causes a larger variance in the realizations of the stochastic process. Next, we will check the generation of trajectories for a segment at \(T=10\) with the same divisions \(N \in \{2^5, 2^7, 2^9\}\).

Finally, we will look at the situation when \(T=100\) with the same \(N\) divisions.

We can see how an increase in the time segment interval allows the stochastic process to increase the achieved values. Higher deviations still show up for a small division of the segment due to the relatively large size of \(\Delta t\). It is also an interesting observation that for the same reason, the process appears to be more stable for a large section division. Let’s additionally generate several realizations on the same segment and with the same division to see the randomness of this stochastic process. Assume \(T=100\) and \(N=2^7\).

The above graph fully captures the random nature of stochastic processes. We can see that, despite the same distribution of the section, each process is able to reach a completely different value at given points. Nevertheless, the increments of each process come from a distribution with the same variance.

8 Stochastic mathematical models

When modeling various types of phenomena, it is often necessary to take into account the effect of randomness, which can be done in a number of ways. A natural extension of the deterministic model of ordinary differential equations leads to stochastic differential equation models by adding a noise component to the controlling equations of the system. This approach assumes that there is some degree of noise in the dynamics of the process itself.

We consider an autonomous \(m\)-dimensional initial differential problem

\[ \frac{dx(t)}{dt} = f(x(t)), \quad x(t_0) = x_0, \]

where \(x = (x_1, ..., x_m) \in \mathbb{R}^m\), \(f:\mathbb{R}^m \rightarrow \mathbb{R}^m\), \((t_0, x_0) \in \mathbb{R}^{m+1}\). If the parameters or coefficients in the deterministic model are subject to some random fluctuations then we can include a noise factor which will lead to a system of stochastic differential equations (\(m\)-dimensional It\(\hat{o}\) equation):

\[ dX(t) = f(X(t))dt + g(X(t))dB(t), \quad X(t_0) = X_0, \]

where \(f\) is the deterministic part i.e. coming from a deterministic system and \(g:\mathbb{R}^m \rightarrow \mathbb{R}^m\) \(dB\) is the stochastic perturbation. Furthermore, \(B\) is the \(m\)-dimensional standard Brownian motion process, and \(dB\) is the derivative of It\(\hat{o}\). The solution of \({X(t), t \in [t_0, T]}\) of a stochastic problem is called an It\(\hat{o}\) process.

We can write the differential form in integral form

\[ X(t) = X(t_0) + \int_{t_0}^{t}{f(X(s))ds} + \int_{t_0}^{t}{g(X(s))dB(s)}. \]

As in the deterministic case, we are not always able to determine the exact solution to the problem. We must then use numerical methods that allow us to approximately observe the dynamics of the process \(X\).

In the construction of numerical methods, we will consider an equidistant discretization of the time interval \([t_0, T]\) ie:

\[ t_n = t_0 + n\Delta t, \quad \Delta t = \frac{T - t_0}{N}, \quad n=0,1, ..., N, \]

where \(N\) is a sufficiently large natural number defining the number of subintervals of a segment \([t_0, T]\).

8.1 Euler-Maruyama approximation method.

The simplest method of approximation of the It\(\hat{o}\) process is the Euler method called the Euler-Maruyama method. The Euler-Maruyama approximation of a process \(X\) is called a continuous stochastic process \(\{Y(t), t\in [t_0, T]\}\) satisfying a recursive equation with a deterministic initial condition i.e.:

\[ Y_{n+1} = Y_n + f(Y_n)\Delta t + g(Y_n)\Delta B_n, \quad n=0,1,...,N-1 \\ Y_0 = x_0, \]

where \(Y_n = Y(t_n)\), \(\Delta B_n = B(t_n) - B(t_n - \Delta t) = \sqrt{\Delta t}\epsilon_n\), \(n=1,...,N\). Between the designated discretization points where we have no information about the process \(Y\), the process It\(\hat{o}\) is approximated by linear interpolation, when \(g\equiv 0\) then the method reduces to a deterministic Euler scheme for the differential equation.

8.2 Milstein’s approximation method

Another method is the Milstein method, which has greater accuracy than the Euler-Maruyama method, since it additionally uses a second-order expression from the It\(\hat{o}\)-Taylor expansion around a point for the stochastic case.

The Milstein approximation of the process \(X\) is called a continuous stochastic process \(\{Z(t), t\in [t_0, T]\}\) satisfying a recursive equation with a deterministic initial condition, viz:

\[ Z_{n+1} = Z_n + f(Z_n)\Delta t + g(Z_n)\Delta B_n + \frac{1}{2}g(Z_n)g_X(Z_n)((\Delta B_n)^2 - \Delta t), \\ Z_0 = x_0, \]

where \(n=0,1,...,N-1\) and \(g_X = \frac{dg}{dX}\).

8.3 Simulation numerical methods

We will consider the process It\(\hat{o}\) \(\{X(t), t\in [0,T]\}\) satisfying the linear stochastic initial problem.

\[ dX(t) = aX(t)dt + bX(t)dB(t), \quad X(0) = X_0 \]

for \(t(t)\in [0,T]\). In this case, \(f=aX(t)\) and \(g(X(t)) = bX(t)\), where \(a\), \(b\) are constants. The solution can be represented in exact form

\[ X(t) = X_0 \exp((a - \frac{1}{2}b^2)t + dB(t)). \]

We will compare the exact solution determined in the discrete points with the form of the trajectory of the process \(Y\) obtained by the Euler method, i.e.:

\[ X(t_n) = X_0 \exp((a-\frac{1}{2}b^2)t_n + b\sum_{i=1}^{n}{\Delta B_{i-1}}). \] To do this, we will transform the function from the previous lab to additionally return a vector of values \(\Delta B_n = \sqrt{\Delta t}\epsilon_n\) for \(n=1, ..., N\).

brown <- function(N, T, s){
  set.seed(s)
  b <- cumsum(c(0, sqrt(T/N) * rnorm(N, 0, 1)))
  b_d <- sqrt(T/N) * rnorm(N, 0, 1)
  return(list(B = b, delta_B = b_d))
}

WLet’s now generate the difference vector along with the trajectory of the Brownian motion process for \(T=1\) and for \(N=2^m\), where \(m=12\) at a fixed \(\omega\), i.e. a single realization of the process.

For the values thus generated and the parameters assumed, we will find the exact solution by assuming the following

\[ (a, b, X_0) = (1.5, 1, 1). \]

The results are shown below in the chart .

We will then generate the trajectory of the Brownian motion process for \(N=2^8\) and present it in a single graph with the exact solution \(X(t_n)\) for different noise parameters \(b\) taken from the interval \([0,1]\) with a step of \(0.1\).

We can see how the noise parameter affects the stochasticity of a given process. Indeed, it regulates the effect of Brownian motion trajectories on the exact solution by which the larger the value of the parameter, the larger the deviation.

We will now generate Brownian motion trajectories for \(N=2^8\) and \(T=1\) and the parameter set as before, i.e. \((a, b, X_0) = (1.5, 1, 1)\).

Let us also introduce an approximation algorithm for the \(X\) Euler-Maruyama process. This approximation is called the stochastic process \(\{Y(t), t\in [t_0, T]\}\) satisfying the recursive equation

\[ \begin{cases} Y_{n+1} = Y_n + f(Y_n)\Delta t + g(Y_n)\Delta B_n, \: n=1, ..., N-1, \\ Y_0 = x_0, \end{cases} \]

gdzie \(Y_n = Y(t_n)\), \(\Delta B_n = B(t_n) - B(t_n - \Delta t) = \sqrt{\Delta t}\epsilon_n\), \(n=1, ..., N\).

Let us also introduce Milstein’s approximation algorithm. Such an approximation is called a stochastic process \(\{Z(t), t\in [t_0, T]\}\) satisfying the recursive equation

\[ \begin{cases} Z_{n+1} = Z_n + f(Z_n)\Delta t + g(Z_n)\Delta B_n + \frac{1}{2}g(Z_n)g_X(Z_n)((\Delta B_n)^2 - \Delta t) \\ Z_0 = x_0 \end{cases} \]

for \(n=0, ..., N-1\), \(g_X = \frac{\partial g}{\partial X}\).

In our case, considering the process It\(\hat{o}\) \(X(t)\) satisfying a linear stochastic differential equation

\[ dX(t) = aX(t) dt + bX(t)dB(t), \: X(0)=x_0 \]

In the denotations from the algorithms, we take \(f(X(t)) = aX(t)\) and \(g(X(t)) = bX(t)\) where \(a\) and \(b\) are constants. For this reason, we will represent the trajectories of the process \(Y(t)\) obtained by the Euler-Maruyama method as an iterative solution of the

\[ Y_{n+1}=Y_n + aY_{n}\Delta t + bY_{n}\Delta B_n \]

and in Milstein’s method the trajectories of the process \(Z\) as an iterative solution

\[ Z_{n+1} = Z_n + aZ_{n}\Delta t + bZ_n \Delta B_n + \frac{1}{2}bZ_n\cdot b\cdot ((\Delta B_n)^2 - \Delta t). \]

Using these algorithms, we will determine the trajectories of the process \(Y_n\), \(Z_n\) for \(N\in \{2^2, 2^4, 2^6, 2^8\}\) and contrast them with the trajectory of \(X_n\).

The trajectories of the exact solution are shown by the bold black line. We can observe how both methods approximate the real solution very accurately at the same pitch. Although by decreasing the value of \(N\) we observe an increasingly poorer fit, the methods still reflect the dynamics of the real solution, and the calculated estimates are surprisingly close to the true values. Let’s repeat the task for another implementation.

As you can see, the situation is analogous. Let’s consider one more realization.

Again, the same conclusions. So indeed the obtained approximations are very accurate at the same pitch. When the value of \(N\) is low then, as mentioned, the approximations are significantly more outliers than the actual solutions but nevertheless shape fits the trajectory of the exact solution and the difference is mainly on the value reached by the process. Let us now present the exact solution and its approximations on the same value of the division. In the first place for \(N=2^5\).

We can observe the same shape of the trajectory with slight deviations of the approximation from the exact solution. Let’s check \(N=2^6\).

Here juz these deviations practically disappear and it is hard to find differences. Only at the end of the interval they begin to reveal themselves. Let’s check \(N=2^{10}\).

And in this situation the exact solution is already almost completely covered by the approximations. The large scale makes it a little difficult to interpret the results, but we can see that while the Milstein method was noticeable in the previous situation, in this one it is completely “covered” by the exact solution. This suggests better accuracy than the Euler-Milstein method, which is due to an additional factor in the form of a partial derivative after the function responsible for the stochastic growth in the process.

All the simulations carried out show the enormity of the possibilities of numerical methods and the potential of this. Solving a stochastic differential equation is not an easy task, and thanks to numerical methods with high correctness we are able to estimate these solutions knowing only the original form of the equation.

9 Deterministic vs stochastic mathematical models

We consider a stochastic predator-prey mathematical model of the form \[ \begin{cases} dP(t) = [rP(t)(1-\frac{P(t)}{K}) - sP(t)Q(t)]dt + \sigma_1P(t)B(t),\\ dQ(t) = [-uQ(t) + asP(t)Q(t)]dt + \sigma_2 Q(t)B(t). \end{cases} \]

We assume that the strength of influencing populations from the environment is the same therefore \(\sigma_1 = \sigma_2 = \sigma\). Note that if \(\sigma =0\) then we get the deterministic mathematical model analyzed in earlier examples. Models based on stochastic differential equations fit the reality of the modeled phenomena better than deterministic models. This is because in reality the phenomena are affected by many random, unpredictable factors, the effect of which can accumulate over time, significantly affecting the dynamics of the phenomenon. This, in turn, will cause huge differences between stochastic and deterministic models. There are situations when the solutions of deterministic models stablize at a certain level, and stochastic models on the contrary - with a high coefficient of randomness are extremely irregular and unpredictable. For this reason, it is extremely important to analyze the phenomenon under study and take into account the appropriate level of random noise.

When modeling stochastically, we must take into account the fact that solutions can “explode”, that is, they can become unbounded over a finite stretch of time. Mostly this happens when the conditions of the existence theorem and unambiguous growth are not satisfied (among other things, the functions contained in the stochastic differential equation should be Lipschitz and satisfy the linear growth condition).

We will now conduct an analysis of the predator-prey model in question for a set of parameters

\[ (P_0, Q_0, r, u, s, a, K) = (0.9, 0.1, 0.5, 0.5, 1.5, 1.1, 1), \quad (T, N) = (80, 10000) \]

for different values of \(\sigma\) to conduct an analysis of the effect of the stochastic factor on the deterministic model.

First, let’s plot the corresponding trajectories by value \(\sigma\).

In the example above, we have shown the effect of the stochastic factor on the deterministic model. We can see how already a small value of \(\sigma\)-noise affects the behavior of the trajectory and introduces irregularities. As the value of noise increases, the deviation from the deterministic solution increases. Already for \(\sigma = 0.05\) the oscillations appear in other places and the cyclicity of the solutions changes in a sense. Nevertheless, in the phase portrait we can still see how the solution strives towards the deterministic solution, but there are large disturbances in its neighborhood. In general, random noise has a big impact on changing the dynamics of both populations. In the deterministic scenario, we have small fluctuations and rapid convergence, while when the stochastic factor is introduced, both populations experience rapid turbulence, and stabilization is achieved much later. The phenomena of rapid population growth, as well as rapid extinction, can be seen in certain cycles.

However, the above results are obtained from a single Brownian motion trajectory. To better observe the long-term behavior of the stochastic model, we will perform simulations for the same parameters but different Brownian motion trajectories. In addition, so that the graphs remain readable we will analyze different values of \(\sigma\) noise on separate graphs. We will start with a value of \(\sigma = 0.01\).

We can observe that at this value of noise, any periodicity disappears, and any Brownian motion generated does not manifest any patterned behavior relative to the solution trajectory. Nevertheless, an approximate convergence to a deterministic solution is achieved, which we can observe on the phase portrait. On the solution graphs, too, it is evident that in the long run, when the stochastic factor is taken into account, we will be close to the equilibrium point of this system, and we could say that it reflects to the behavior of the deterministic model. Still, as in the example above, there are small oscillations, and the dynamics of the two populations are not clearly defined. We can see how the abundance of both predators and prey increases rapidly in some places, only to then decrease just as rapidly. Let’s raise the noise value to \(\sigma = 0.05\).

The solutions have changed drastically despite the relatively small increase in noise values. The solution trajectories show significantly increased amplitudes with increased irregularity compared to the earlier example. In the prey population, these are noticeable to a greater extent than in the predator population. In the phase portrait, we can see how solutions with a stochastic factor begin to diverge from the deterministic solution. Still, with much larger fluctuations, they are eventually around a certain equilibrium point considering the higher level of deviation. Nevertheless, we are no longer able to say that these solutions reflect the behavior of the deterministic model. We can still see periodic increases in abundance occurring alternately with extinction events among populations. Let’s increase the noise factor more dramatically to a level of \(\sigma=0.5\).

Here we are no longer able to see the behavior of the deterministic model to any degree. The amplitudes of the solutions have increased dramatically without the slightest regularity. Moreover, the solutions do not seem to disappear over time as they did in the previous examples. A much higher frequency of perturbations appeared in the predator population than in the prey population where huge fluctuations are less frequent. In the phase portrait, we see no convergence to a certain equilibrium point, as the values obtained have become completely random. The above example illustrates how an excessive degree of randomness in the model can completely muddy the actual deterministic solution.

In the next step, we will examine the distribution of \(1000\) sample solutions of the stochastic model at the end of the time interval. We will consider numerical parameters \((T, N) = (50, 1000)\) for this purpose, so we will generate \(1000\) trajectories of solutions for the values of the noise parameter \(\sigma \in \{0.05, 0.1, 0.3, 0.5\}\) and the model parameters from the previous examples.

We can see that for low values of the noise \(\sigma \in \{0.05, 0.1\}\) the values accumulate around the point \(P^* \approx 0.3\) for \(P\) solutions and \(Q^* \approx 0.25\) for \(Q^*\) solutions. If we go back to previous analyses, including the first visualization of the trajectory, we see that these points are approximations of the equilibrium points of the deterministic model under consideration. Thus, the solution will still converge to the equilibrium point of the deterministic model in most cases, so in a sense it replicates its behavior. However, we cannot say this for higher values of \(\sigma \in \{0.3, 0.5\}\). Then, for \(P\) solutions, the points are more spread out over the entire interval \([0, 0.5]\), and for \(Q\) solutions they are concentrated at the zero value. Such a scenario would mean that the predator population would remain at a certain level and the prey population would die out, which is not consistent with the interpretation of the Predator-Prey model, since predators need food to survive. So here, too, the excessive randomness factor completely disturbs the model’s behavior and the total distribution interpretation is consistent with the earlier analysis.

Bibliography

1.

Szafrańska A (2022) Mathematical Modeling and Computer Simulations - Lectures.

2.

Strogatz SH (1988) Love affairs and differential equations.

3.

Jury EI (1971) ‘Inners’approach to some problems of system theory.